1 Anàlisi de dissenys factorials de dos o més factors creuats

Podeu accedir a una versió electrònica actualitzada d’aquest document clicant al següent enllaç aminarroa.github.io/DEAD.html

1.1 Anova de dos factors fixos

1.1.1 Exemple 1: Analitzant dues causes de variabilitat

mouse

Un investigador vol avaluar l’eficàcia de 4 possibles vacunes per a una malaltia. L’investigador opina també que el número de dosis subministrades: una, dos o tres, pot tenir una influència en el nivell de protecció aconseguit. Efectuant un estudi en rates de laboratori opta per dividir un total de 48 rates en grups de 4 i sotmetre a cada grup a una vacuna i un número de dosis prefixades. Es mesura la quantitat d’anticòs Ig G en gr/l present a l’organisme.

anticos<-c(6.1,8.2,6.3,6.5,7.5,11.0,9.5,7.1,6.6,8.8,6.3,6.6,4.3,5.2,5.6,6.2,
           3.6,5.2,4.4,5.6,2.9,5.1,3.5,10.2,4.0,4.9,3.1,7.1,2.3,12.4,4.0,3.8,
           2.2,3.0,2.3,3.0,2.1,3.7,2.5,3.6,1.8,3.8,2.4,3.1,2.3,2.9,2.2,3.3)
vacuna<-gl(4,12)
dosis<-factor(rep(1:3,times=16))
immuno<-data.frame(anticos,vacuna,dosis)
  • L’objectiu és estudiar simultàniament si els dos factors vacuna i dosi, amb 4 i 3 nivells respectivament, poden influir en la variable resposta quantitativa quantitat d’anticòs.

  • Els individus han estat aleatoritzats respecte els factors.

1.1.2 Disseny de dos factors amb interacció. Estructura de les dades

L’estructura de les dades correspon a una taula bidimensional on cada cel·la equival a una combinació de nivells dels dos factors, i dintre de cada cel·la tenim les rèpliques que corresponen a la condició experimental que representa la combinació.

\[ \begin{array}{cccc} & B_1 & \cdots & B_b \\ A_1 & y_{111},...,y_{11n_{11}} & \cdots & y_{1b1},...,y_{1bn_{1b}}\\ \vdots & \vdots & & \vdots \\ A_a & y_{a11},...,y_{a1n_{a1}} & \cdots & y_{ab1},...,y_{abn_{ab}} \end{array} \]

Si el disseny és balancejat \(n_{ij} = n\) per tot \(i=1,...,a\) i \(j=1,...,b\).

D’ara endavant si no diem el contrari assumirem que el disseny és balancejat.

1.1.3 Model lineal de l’ANOVA de dos factors fixos creuats

El model lineal que respon a l’enunciat de l’experiment és:

\[ Y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_{ij} +\epsilon_{ijk} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,n_{ij} \]

on

  • \(Y_{ijk}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\),

  • \(\mu\) és la mitjana general,

  • \(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor A,

  • \(\beta_j\) és un paràmetre del model que mesura l’efecte del nivell \(j\) del factor B,

  • \(\gamma_{ij}\) és un paràmetre del model que mesura l’efecte de la interacció del nivell \(i\) del factor A amb el nivell \(j\) del factor B (també és freqüent representar-lo com \((\alpha \beta)_{ij})\),

  • \(\epsilon_{ijk}\) és l’error aleatori corresponent a l’observació \(i,j,k\) i dels quals suposem que són independents i idènticament distribuïts amb \(E(\epsilon_{ijk}) = 0\) i \(Var(\epsilon_{ijk})=\sigma^2\).

  • Assumirem també les següents restriccions sobre els paràmetres \(\sum_{i=1}^a \alpha_i = 0\), \(\sum_{j=1}^b \beta_j = 0\), \(\sum_{i=1}^a \gamma_{ij} = 0 \,\, j=1,...,b\) i finalment \(\sum_{j=1}^b \gamma_{ij} = 0 \,\, i=1,...,a\).

Com veurem serà necessari també assumir una distribució normal sobre els errors

\[ \epsilon_{ijk} \sim N(0,\sigma)\]

i per tant cada observació experimental

\[ Y_{ijk} \sim N(\mu_{ij},\sigma) \]

on \(\mu_{ij} = \mu + \alpha_i + \beta_j + \gamma_{ij}\).

La interacció la podem interpretar com l’efecte de la combinació de dos nivells, un de cada factor, que provoca una modificació en la variable resposta que va més enllà de l’efecte merament additiu.

1.1.4 Contrast d’hipòtesis

L’Anova de dos factors creuats amb interacció permet contrastar les següents hipòtesis:

  • Hipòtesi sobre la significació del factor A

\[ H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0 \] \[ H_1: \alpha_i \ne \alpha_{i'} \mbox{ per algun } i \ne i' \,\,\, ( \mbox{ equivalentment }H_1: \exists \alpha_i \ne 0)\]

  • Hipòtesi sobre la significació del factor B

\[ H_0: \beta_1 = \beta_2 = \cdots = \beta_b = 0 \] \[ H_1: \beta_j \ne \beta_{j'} \mbox{ per algun } j \ne j' \,\,\, ( \mbox{ equivalentment } H_1: \exists \beta_j \ne 0)\]

  • Hipòtesi sobre la significació de la interacció

\[ H_0: \gamma_{11} = \gamma_{12} = \cdots = \gamma_{ab} = 0 \] \[ H_1: \exists \gamma_{ij} \ne 0\]

1.1.5 Estimacions dels paràmetres del model

1.1.5.1 Definicions

\(N=abn \,\,\,\) (Disseny balancejat)

\(Y_{i..} = \sum_{j=1}^{b} \sum_{k=1}^{n}Y_{ijk}\)

\(Y_{.j.} = \sum_{i=1}^{a} \sum_{k=1}^{n}Y_{ijk}\)

\(Y_{ij.} = \sum_{k=1}^{n} Y_{ijk}\)

\(Y_{...} = \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} Y_{ijk}\)

\(\overline{Y}_{i..}= \frac{Y_{i..}}{bn}\)

\(\overline{Y}_{.j.}= \frac{Y_{.j.}}{an}\)

\(\overline{Y}_{ij.}= \frac{Y_{ij.}}{n}\)

\(\overline{Y}_{...}= \frac{Y_{...}}{N}\)

1.1.5.2 Estimadors puntuals dels paràmetres del model

\(\hat{\mu} = \overline{Y}_{...}\)

\(\hat{\alpha}_i = \overline{Y}_{i..} - \overline{Y}_{...} \,\,\,\, (\alpha_i = \mu_i - \mu)\)

\(\hat{\beta}_j = \overline{Y}_{.j.} - \overline{Y}_{...} \,\,\,\, (\beta_j = \mu_j - \mu)\)

\(\hat{\gamma}_{ij} = \overline{Y}_{ij.} - \overline{Y}_{i..}- \overline{Y}_{.j.}+ \overline{Y}_{...} \,\,\,\, (\gamma_{ij} = \mu_{ij} -\mu_i -\mu_j + \mu)\)

1.1.5.3 Valor predit i residu

\(\hat{Y}_{ijk} = \overline{Y_{ij.}}\)

\(\hat{e}_{ijk} = Y_{ijk} - \overline{Y_{ij.}}\)

1.1.6 Descomposició de la variabilitat

Mantenim la idea de la descomposició de la variabilitat total de les dades. En aquest cas la variabilitat total es descompon en: una part que representa la variabilitat deguda al factor A, una part que correspon a la variabilitat deguda al factor B, una part que correspon a la variabilitat deguda a la interacció dels factors A i B i finalment una part que és la variabilitat dintre dels tractaments o variabilitat del error o residual.

Posteriorment decidirem si els possibles efectes són significatius comparant les variabilitats dels efectes amb la variabilitat residual.

\[ \begin{aligned} SS_T = \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} (Y_{ijk}-\overline{Y}_{...} )^2 = \\ bn \sum_{i=1}^{a} (\overline{Y}_{i..}-\overline{Y}_{...} )^2 + an \sum_{j=1}^{b} (\overline{Y}_{.j.}-\overline{Y}_{...} )^2 + \\ n \sum_{i=1}^{a} \sum_{j=1}^{b} (\overline{Y}_{ij.}-\overline{Y}_{i..}-\overline{Y}_{.j.}+\overline{Y}_{...} )^2 + \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} (Y_{ijk}-\overline{Y}_{ij.} )^2 = \\ SS_A + SS_B + SS_{AB} + SS_E \\ \end{aligned} \]

on

  • \(SS_T\) recull la variabilitat total de les dades,

  • \(SS_A\) recull la variabilitat deguda al factor A,

  • \(SS_B\) recull la variabilitat deguda al factor B,

  • \(SS_{AB}\) recull la variabilitat deguda a la interacció,

  • \(SS_E\) recull la variabilitat residual.

1.1.7 Esperances dels quadrats mitjans

1.1.7.1 Quadrats mitjans del error

Si anomenem \(MS_E = \frac{SS_E}{ab(n-1)}\) (quadrats mitjans de l’error) pot demostrar-se que sempre \(E(MS_E) = \sigma^2\).

1.1.7.2 Quadrats mitjans deguts al tractament A

Per una altra banda si anomenem \(MS_A = \frac{SS_A}{a-1}\) (quadrats mitjans del tractament A) pot demostrar-se que \(E(MS_A) = \sigma^2 + \frac{ bn \sum_{i=1}^a \alpha_i^2}{a-1}\).

Per tant i si la hipòtesi nul·la referent al factor \(A\) és certa \((H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0) \Rightarrow E(MS_A)=E(MS_E)=\sigma^2\) i per tant si plantegem

\[ F = \frac{MS_A}{MS_E} \approx 1. \]

En canvi si \(H_0\) és falsa la tendència ha de ser que \(F > 1\).

\[ F = \frac{MS_A}{MS_E} = \frac{\frac{SS_A}{a-1}}{\frac{SS_E}{ab(n-1)}}. \]

  • Si \(H_0\) és certa , \(F\) tendirà a ser propera a 1.

  • Si \(H_0\) és falsa , \(F\) tendirà a ser més gran que 1.

Per tant un criteri raonable per decidir consisteix en rebutjar \(H_0\) si \(F\) és prou gran.

Si \(F \ge c_\alpha \Rightarrow\) rebutgem \(H_0\), on \(c_\alpha\) és triat de forma que \(P_{H_0} (F \ge c_\alpha) = \alpha\).

Per trobar \(c_\alpha\) pot demostrar-se que sota la hipòtesi de errors normals i homocedàstics

\[ \epsilon_{ijk} \sim N(0,\sigma).\]

Sota la hipòtesi nul·la \(MS_A\) i \(MS_E\) segueixen una distribució \(\chi^2\) amb \(a-1\) i \(ab(n-1)\) graus de llibertat respectivament i són independents, i com a conseqüència

\[ F = \frac{MS_A}{MS_E} = \frac{\frac{SS_A}{a-1}}{\frac{SS_E}{ab(n-1)}} \sim F(a-1,ab(n-1))\]

on \(F(a-1,ab(n-1))\) és una distribució \(F\) de Fisher amb \(a-1\) graus de llibertat en el numerador i \(ab(n-1)\) graus de llibertat en el denominador.

1.1.7.3 Quadrats mitjans deguts al tractament B

Anàlogament si anomenem \(MS_B = \frac{SS_B}{b-1}\) pot demostrar-se que \(E(MS_B) = \sigma^2 + \frac{ an \sum_{j=1}^b \beta_j^2}{b-1}\).

Per tant i si la hipòtesi nul·la referent al factor \(B\) és certa \((H_0: \beta_1 = \beta_2 = \cdots = \beta_b = 0) \Rightarrow E(MS_B)=E(MS_E)=\sigma^2\) i per tant si plantegem

\[ F = \frac{MS_B}{MS_E} \approx 1 .\]

En canvi si \(H_0\) és falsa la tendència ha de ser que \(F > 1\).

Sota la hipòtesi nul·la \[ F = \frac{MS_B}{MS_E} = \frac{\frac{SS_B}{b-1}}{\frac{SS_E}{ab(n-1)}} \sim F(b-1,ab(n-1)).\]

1.1.7.4 Quadrats mitjans deguts a la interacció

Si anomenem \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) pot demostrar-se que \(E(MS_{AB}) = \sigma^2 + \frac{ n \sum_{i=1}^a \sum_{j=1}^b \gamma_{ij}^2}{(a-1)(b-1)}\).

Per tant i si la hipòtesi nul·la referent a la interacció és certa \(H_0: \gamma_{11} = \gamma_{12} = \cdots = \gamma_{ab} = 0 \Rightarrow E(MS_{AB})=E(MS_E)=\sigma^2\) i per tant si plantegem

\[ F = \frac{MS_{AB}}{MS_E} \approx 1. \]

En canvi si \(H_0\) és falsa la tendència ha de ser que \(F > 1\).

Sota la hipòtesi nul·la \[ F = \frac{MS_{AB}}{MS_E} = \frac{\frac{SS_{AB}}{(a-1)(b-1)}}{\frac{SS_E}{ab(n-1)}} \sim F((a-1)(b-1),ab(n-1)).\]

1.1.8 Taula de l’Anàlisi de la variància de dos factors (ANOVA)

Font de variació Suma de quadrats Graus de llibertat Quadrats mitjans Estadístic F
Factor A \(SS_A\) \(a-1\) \(MS_A = \frac{SS_A}{a-1}\) \(F = \frac{MS_A}{MS_E}\)
Factor B \(SS_B\) \(b-1\) \(MS_B = \frac{SS_B}{b-1}\) \(F = \frac{MS_B}{MS_E}\)
Interacció AB \(SS_{AB}\) \((a-1)(b-1)\) \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) \(F = \frac{MS_{AB}}{MS_E}\)
Error \(SS_E\) \(ab(n-1)\) \(MS_E = \frac{SS_E}{ab(n-1)}\)
Total \(SS_T\) \(N-1\)
  • Criteri de decisió

– Calcular el valor dels estadístics \(F\) de la taula a partir de les dades.

– Determinar els valors crítics \(F_{\alpha} (gl_{num},gl_{den})\) per a les corresponents distribucions \(F\) de Fisher d’acord amb \(\alpha\).

Rebutjarem cada hipòtesi nul·la si l’estadístic \(F\) respectiu és més gran que el valor crític \(F \ge F_{\alpha} (gl_{num},gl_{den})\) corresponent o, equivalentment si obtenim el resultat amb un software estadístic, rebutjarem si el p-valor és menor que el nivell de significació.

1.1.9 Resolució amb R

La sintaxi vacuna*dosis representa la generació dels factors principals i la interacció.

## anova
immuno.aov<-aov(anticos~vacuna*dosis,immuno)
anova(immuno.aov)
Df Sum Sq Mean Sq F value Pr(>F)
vacuna 3 140.557292 46.8524306 11.3339325 0.0000221
dosis 2 7.937917 3.9689583 0.9601189 0.3924335
vacuna:dosis 6 4.362083 0.7270139 0.1758698 0.9816806
Residuals 36 148.817500 4.1338194 NA NA

En aquest cas direm que hi han diferències degudes a les vacunes, no a les dosis i tampoc hi ha una interacció significativa.

En forma alternativa, sense utilitzar l’asterisc i especificant tots els termes:

immuno.aov<-aov(anticos~vacuna+dosis+vacuna:dosis,immuno)

Gràfica descriptiva: Interaction plot

interaction.plot(immuno$vacuna,immuno$dosis,immuno$anticos,col='red',trace.label='dosis',xlab='Vacuna',ylab='Anticos')

1.1.9.1 Taules de coeficients i mitjanes

## Estimació de mitjanes
taula<-model.tables(immuno.aov,type='mean')
taula
## Tables of means
## Grand mean
##          
## 4.960417 
## 
##  vacuna 
## vacuna
##     1     2     3     4 
## 7.542 5.150 4.342 2.808 
## 
##  dosis 
## dosis
##     1     2     3 
## 5.000 4.444 5.437 
## 
##  vacuna:dosis 
##       dosis
## vacuna 1     2     3    
##      1 7.725 7.275 7.625
##      2 5.000 4.475 5.975
##      3 4.525 3.325 5.175
##      4 2.750 2.700 2.975
## Estimació d'efectes (paràmetres del model)
model.tables(immuno.aov)
## Tables of effects
## 
##  vacuna 
## vacuna
##       1       2       3       4 
##  2.5812  0.1896 -0.6187 -2.1521 
## 
##  dosis 
## dosis
##       1       2       3 
##  0.0396 -0.5167  0.4771 
## 
##  vacuna:dosis 
##       dosis
## vacuna 1       2       3      
##      1  0.1438  0.2500 -0.3938
##      2 -0.1896 -0.1583  0.3479
##      3  0.1437 -0.5000  0.3563
##      4 -0.0979  0.4083 -0.3104
  • Algunes estimacions dels paràmetres del model:

\(\hat{\mu} =\) 4.9604167

\(\hat{\alpha}_1\) = 7.5416667 - 4.9604167= 2.58125

\(\hat{\beta}_2\) = 4.44375 - 4.9604167= -0.5166667

\(\hat{\gamma}_{23}\) = 5.975 - 5.15 - 5.4375 + 4.9604167= 0.3479167

\(\hat{\sigma}^2\) = 4.1338194 \(\Rightarrow \hat{\sigma}\) = 2.0331796

1.1.10 Comparacions múltiples

Assumint que algun dels factors principals ha estat significatiu i la interacció no, podem plantejar-nos les comparacions múltiples per a detectar quins nivells del factor principal són diferents entre ells.

library(agricolae)
HSD.test(immuno.aov,'vacuna',console=T)
## 
## Study: immuno.aov ~ "vacuna"
## 
## HSD Test for anticos 
## 
## Mean Square Error:  4.133819 
## 
## vacuna,  means
## 
##    anticos       std  r        se Min  Max   Q25  Q50   Q75
## 1 7.541667 1.5382743 12 0.5869284 6.1 11.0 6.450 6.85 8.350
## 2 5.150000 1.8705857 12 0.5869284 2.9 10.2 4.125 5.15 5.600
## 3 4.341667 2.8833246 12 0.5869284 2.2 12.4 2.825 3.45 4.225
## 4 2.808333 0.6841828 12 0.5869284 1.8  3.8 2.275 2.70 3.375
## 
## Alpha: 0.05 ; DF Error: 36 
## Critical Value of Studentized Range: 3.808798 
## 
## Minimun Significant Difference: 2.235492 
## 
## Treatments with the same letter are not significantly different.
## 
##    anticos groups
## 1 7.541667      a
## 2 5.150000      b
## 3 4.341667     bc
## 4 2.808333      c

En aquest cas no seria correcte fer les comparacions múltiples pel factor dosis perquè no ha estat significatiu.

Però, com a exercici didàctic, què passa si ho fem:

HSD.test(immuno.aov,'dosis',console=T)
## 
## Study: immuno.aov ~ "dosis"
## 
## HSD Test for anticos 
## 
## Mean Square Error:  4.133819 
## 
## dosis,  means
## 
##   anticos      std  r        se Min  Max   Q25  Q50   Q75
## 1 5.00000 2.196968 16 0.5082949 2.1  9.5 3.450 4.35 6.275
## 2 4.44375 2.027468 16 0.5082949 1.8  8.2 2.900 3.75 5.775
## 3 5.43750 3.262693 16 0.5082949 2.2 12.4 2.975 4.50 6.600
## 
## Alpha: 0.05 ; DF Error: 36 
## Critical Value of Studentized Range: 3.456758 
## 
## Minimun Significant Difference: 1.757053 
## 
## Treatments with the same letter are not significantly different.
## 
##   anticos groups
## 3 5.43750      a
## 1 5.00000      a
## 2 4.44375      a

Tots els resultats ens permeten reinterpretar el gràfic d’interaccions.

interaction.plot(immuno$vacuna,immuno$dosis,immuno$anticos,col='red',trace.label = 'dosis')

1.1.11 Interpretació de la interacció

A l’exemple anterior hem vist que la interacció no era significativa, això es traduïa en un paral·lelisme en el gràfic d’interaccions

interaction.plot(immuno$vacuna,immuno$dosis,immuno$anticos,col='red',trace.label = 'dosis')

Al següent apartat veurem un altre exemple força diferent.

1.1.11.1 Exemple 2: Grau de somnolència en fàrmacs analgèsics

En fàrmacs per reduir el dolor, el grau de somnolència és un dels efectes secundaris que de forma habitual s’avalua en base a la consideració del temps de reacció a un estímul (en segons). En aquesta direcció, estem interessats en contrastar l’efecte de tres medicaments (A, B i C) i en la forma d’administrar-los (C: càpsula, I: injectable). Els resultats obtinguts en individus de característiques similars són els següents

temps<-c(4.35,3.86,4.62,5.14,4.8,3.71,3.9,3.4,5.97,6.05,6.32,
         6.07,1.35,1.62,0.34,0.63,1.08,1.51,1.61,0.94,0.84,0.81,1.67,0.96)
admin<-gl(2,12)
farmac<-factor(rep(rep(1:3,each=4),2))
somnol<-data.frame(temps,admin,farmac)
interaction.plot(somnol$farmac,somnol$admin,somnol$temps,col='red',trace.label = 'admin')

somnol.aov<-aov(temps~farmac*admin,somnol)
anova(somnol.aov)
Df Sum Sq Mean Sq F value Pr(>F)
farmac 2 4.449900 2.2249500 10.27420 0.0010553
admin 1 83.738704 83.7387042 386.68215 0.0000000
farmac:admin 2 5.749633 2.8748167 13.27511 0.0002869
Residuals 18 3.898025 0.2165569 NA NA

Quan la interacció és significativa això ens vol dir que els efectes principals no són independents un de l’altre i per tant es fa difícil d’interpretar el sentit de les possibles diferències entre els nivells dels factors. Es més recomanable analitzar les diferències entre els nivells de la interacció o fer comparacions per a un factor dins cada nivell de l’altre factor.

HSD.test(somnol.aov,c('admin','farmac'),console=T)
## 
## Study: somnol.aov ~ c("admin", "farmac")
## 
## HSD Test for temps 
## 
## Mean Square Error:  0.2165569 
## 
## admin:farmac,  means
## 
##      temps       std r        se  Min  Max    Q25   Q50    Q75
## 1:1 4.4925 0.5341270 4 0.2326784 3.86 5.14 4.2275 4.485 4.7500
## 1:2 3.9525 0.6014081 4 0.2326784 3.40 4.80 3.6325 3.805 4.1250
## 1:3 6.1025 0.1512999 4 0.2326784 5.97 6.32 6.0300 6.060 6.1325
## 2:1 0.9850 0.5995832 4 0.2326784 0.34 1.62 0.5575 0.990 1.4175
## 2:2 1.2850 0.3252179 4 0.2326784 0.94 1.61 1.0450 1.295 1.5350
## 2:3 1.0700 0.4052160 4 0.2326784 0.81 1.67 0.8325 0.900 1.1375
## 
## Alpha: 0.05 ; DF Error: 18 
## Critical Value of Studentized Range: 4.49442 
## 
## Minimun Significant Difference: 1.045754 
## 
## Treatments with the same letter are not significantly different.
## 
##      temps groups
## 1:3 6.1025      a
## 1:1 4.4925      b
## 1:2 3.9525      b
## 2:2 1.2850      c
## 2:3 1.0700      c
## 2:1 0.9850      c

Podem veure diverses situacions a la següent figura extreta del llibre de G. Quinn i M. Keough Experimental Design and Data Analysis for Biologists, Cambridge University Press, 2002.

1.1.12 Verificació de les premisses del ANOVA

1.1.12.1 Exemple 2: Grau de somnolència en fàrmacs analgèsics

1.1.12.1.1 Normalitat dels residus
shapiro.test(somnol.aov$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  somnol.aov$residuals
## W = 0.95564, p-value = 0.3572
plot(somnol.aov,which=2)

1.1.12.1.2 Homogeneïtat de variàncies
bartlett.test(temps~interaction(admin,farmac),somnol)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  temps by interaction(admin, farmac)
## Bartlett's K-squared = 5.1654, df = 5, p-value = 0.396
plot(somnol.aov,which=1)

1.1.12.2 Exemple 1: vacunes i dosis

1.1.12.2.1 Normalitat dels residus
shapiro.test(immuno.aov$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  immuno.aov$residuals
## W = 0.87414, p-value = 0.0001017
plot(immuno.aov,which=2)

1.1.12.2.2 Homogeneïtat de variàncies
bartlett.test(anticos~interaction(vacuna,dosis),immuno)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  anticos by interaction(vacuna, dosis)
## Bartlett's K-squared = 26.217, df = 11, p-value = 0.006025
plot(immuno.aov,which=1)

En aquest cas veiem clarament tant una desviació de la normalitat com una heteroscedasticitat (variàncies desiguals). Una possible forma de corregir el problema pot ser una transformació de la variable resposta.

Algunes de les transformacions més usuals poden trobar-se a la taula següent:

Tipus de dades Transformació Instrucció R
Mesures continues ln(x) log(x)
- ln(x+1) log(x+1)
- log(x) log(x,10)
Comptatges \(\sqrt{x}\) sqrt(x)
Percentatges arcsin asin(sqrt(x))*180/pi
1.1.12.2.3 Anàlisi amb les dades transformades

Nosaltres hem optat per aplicar la transformació logarítmica en base \(e\).

immuno.aov.trans<-aov(log(anticos)~vacuna*dosis,immuno)
anova(immuno.aov.trans)
Df Sum Sq Mean Sq F value Pr(>F)
vacuna 3 6.4060309 2.1353436 16.3520779 0.0000007
dosis 2 0.2059845 0.1029923 0.7886962 0.4621314
vacuna:dosis 6 0.1281582 0.0213597 0.1635688 0.9848036
Residuals 36 4.7010766 0.1305855 NA NA
shapiro.test(immuno.aov.trans$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  immuno.aov.trans$residuals
## W = 0.94396, p-value = 0.0231
plot(immuno.aov.trans,which=2)

bartlett.test(log(anticos)~interaction(vacuna,dosis),immuno)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  log(anticos) by interaction(vacuna, dosis)
## Bartlett's K-squared = 15.893, df = 11, p-value = 0.1451
plot(immuno.aov.trans,which=1)

Com veiem els problemes de no normalitat i hetersocedasticitat s’han corregit.

1.1.13 Cas d’una rèplica per nivell de la interacció

  • Si tenim una única rèplica per cada combinació de nivells dels dos factors (una rèplica per casella) \(n=1\) la suma de quadrats del error \(SS_E\) té 0 graus de llibertat i \(\sigma^2\) no és estimable.

  • La solució es no incloure la interacció en el model. L’equació del model passa a ser:

\[ Y_{ij} = \mu + \alpha_i + \beta_j +\epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b , \]

i en aquest cas \(SS_T = SS_A + SS_B + SS_E\), adaptant les fórmules al cas \(n=1\).

Formalment la interacció no pot separar-se dels residus.

La taula ANOVA queda

Font de variació Suma de quadrats Graus de llibertat Quadrats mitjans Estadístic F
Factor A \(SS_A\) \(a-1\) \(MS_A = \frac{SS_A}{a-1}\) \(F = \frac{MS_A}{MS_E}\)
Factor B \(SS_B\) \(b-1\) \(MS_B = \frac{SS_B}{b-1}\) \(F = \frac{MS_B}{MS_E}\)
Error \(SS_E\) \((a-1)(b-1)\) \(MS_E = \frac{SS_E}{(a-1)(b-1)}\)
Total \(SS_T\) \(N-1\)

Per verificar en el cas d’una rèplica si és assumible la suposició de no interacció existeixen alguns tests com per exemple el Test de additivitat de Tukey.

1.1.13.1 Subexemple de l’exemple 1: vacunes i dosis

anticos<-c(6.1,8.2,6.3,4.3,5.2,5.6,4.0,4.9,3.1,2.1,3.7,2.5)
vacuna<-gl(4,3)
dosis<-factor(rep(1:3,times=4))
immuno1<-data.frame(anticos,vacuna,dosis)
## anova amb la interacció
immuno1.aov<-aov(anticos~vacuna*dosis,immuno1)
anova(immuno1.aov)
Df Sum Sq Mean Sq F value Pr(>F)
vacuna 3 27.086667 9.0288889 NaN NaN
dosis 2 4.291667 2.1458333 NaN NaN
vacuna:dosis 6 2.288333 0.3813889 NaN NaN
Residuals 0 0.000000 NaN NA NA
## anova sense la interacció
immuno1.aov<-aov(anticos~vacuna+dosis,immuno1)
anova(immuno1.aov)
Df Sum Sq Mean Sq F value Pr(>F)
vacuna 3 27.086667 9.0288889 23.673707 0.0010034
dosis 2 4.291667 2.1458333 5.626366 0.0420611
Residuals 6 2.288333 0.3813889 NA NA

1.2 Dos factors aleatoris

1.2.1 Exemple 3. Variabilitat de la mostra i del analista.

ecoli

Es vol avaluar la variabilitat en la determinació de la qualitat d’un medi de cultiu per tal d’estudiar la rapidesa d’aparició de noves colònies d’Escherichia coli, en base al recompte d’individus en un cultiu durant un període de 48 hores. Es creu que els errors en el comptatge es poden atribuir a la preparació de la mostra o al comptatge per part dels analistes. Per aquest motiu, s’han realitzat 3 preparacions triades a l’atzar de la mateixa base i es realitza l’anàlisi de cada una d’elles per part de dos operaris triats també a l’atzar. Volem valorar la variabilitat aportada a la variable resposta pels diferents factors.

Els resultats obtinguts són:

compt<-c(343,344,341,339,365,364,359,361,321,325,317,322)
prep<-gl(3,4)
analist<-factor(rep(rep(1:2,each=2),3))
dades<-data.frame(prep,analist,compt)
dades
prep analist compt
1 1 343
1 1 344
1 2 341
1 2 339
2 1 365
2 1 364
2 2 359
2 2 361
3 1 321
3 1 325
3 2 317
3 2 322
  • L’objectiu es veure si l’analista o la preparació de la mostra tenen influència en la variable resposta.

  • No estem interessants especialment ni en els dos analistes en concret que han fet el recompte ni en les tres preparacions analitzades, voldríem que les conclusions fossin vàlides pel conjunt de preparacions i analistes possibles, els dos factors són aleatoris.

1.2.2 Model lineal del ANOVA de dos factors aleatoris

Suposem ara que tenim un experiment amb dues fonts de variabilitat controlades \(A\) i \(B\) però que les dues són aleatòries, es a dir els seus nivells són mostres aleatòries de mida \(a\) i \(b\) respectivament d’un conjunt més gran de possibles nivells.

El model lineal és ara

\[ Y_{ijk} = \mu + A_i + B_j + AB_{ij} +\epsilon_{ijk} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,n_{ij} \]

on

  • \(Y_{ijk}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\).

  • \(\mu\) és la mitjana general.

  • \(A_i\) és una variable aleatòria, no un paràmetre constant. Assumim que \(A_i \sim N(0,\sigma_A)\).

  • \(B_j\) és una variable aleatòria, no un paràmetre constant. Assumim que \(B_j \sim N(0,\sigma_B)\).

  • \(AB_{ij}\) és una variable aleatòria, no un paràmetre constant. Assumim que \(AB_{ij} \sim N(0,\sigma_{AB})\).

  • \(\epsilon_{ijk}\) és l’error aleatori corresponent a l’observació \(i,j,k\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ijk} \sim N(0,\sigma)\).

  • Assumirem també que les v.a. \(A_i\), \(B_j\), \(AB_{ij}\) i \(\epsilon_{ijk}\) són v.a. independents.

i per tant cada observació experimental

\[ Y_{ijk} \sim N(\mu,\sqrt{\sigma^2_A+\sigma^2_B+\sigma^2_{AB}+\sigma^2} ).\]

1.2.2.1 Dependència entre les observacions

Tal i com passa en el model aleatori d’un factor, l’existència de factors aleatoris introdueix una dependència entre les observacions dintre d’algunes cel·les del disseny. Pot demostrar-se que:

\[ Cov(Y_{ijk},Y_{i'j'k'}) = \left\{ \begin{array}{ll} \sigma^2_A+\sigma^2_B+\sigma^2_{AB} & \mbox{ si } i=i', j=j', k \ne k' \\ \sigma^2_A & \mbox{ si } i=i', j \ne j' \\ \sigma^2_B & \mbox{ si } i \ne i', j = j' \\ 0 & \mbox{ si } i \ne i', j \ne j' \end{array} \right. \]

Evidentment si \(i=i', j=j', k = k'\) llavors

\[ Cov(Y_{ijk},Y_{ijk})=Var(Y_{ijk})=\sigma^2_A+\sigma^2_B+\sigma^2_{AB}+\sigma^2\]

1.2.3 Contrasts d’hipòtesis

Els contrasts de més interès ara són els següents

  • Variabilitat aportada pel factor A \[ H_0: \sigma_A^2 = 0 \] \[ H_1: \sigma_A^2 > 0\]

  • Variabilitat aportada pel factor B \[ H_0: \sigma_B^2 = 0 \] \[ H_1: \sigma_B^2 > 0\]

  • Variabilitat aportada per la interacció \[ H_0: \sigma_{AB}^2 = 0 \] \[ H_1: \sigma_{AB}^2 > 0\]

La descomposició de la suma de quadrats continua sent vàlida \(SS_A + SS_B +SS_{AB} + SS_E\) i les sumes de quadrats es calculen de la mateixa forma, però els quadrats mitjans tenen un altre sentit i les esperances són ara les següents

  • \(E(MS_E) = \sigma^2\)
  • \(E(MS_A) = \sigma^2 + bn \sigma^2_A+ n \sigma^2_{AB}\)
  • \(E(MS_B) = \sigma^2 + an \sigma^2_B+ n \sigma^2_{AB}\)
  • \(E(MS_{AB}) = \sigma^2 + n \sigma^2_{AB}\)

Per tant els estadístics F per decidir en cadascú dels contrasts anteriors són

\[ F = \frac{MS_A}{MS_{AB}} \sim F(a-1,(a-1)(b-1)) \]

\[ F = \frac{MS_B}{MS_{AB}} \sim F(b-1,(a-1)(b-1)) \]

\[ F = \frac{MS_{AB}}{MS_E} \sim F((a-1)(b-1),ab(n-1)) \]

i la taula ANOVA és la següent

Font de variació Suma de quadrats Graus de llibertat Quadrats mitjans Estadístic F
Factor A \(SS_A\) \(a-1\) \(MS_A = \frac{SS_A}{a-1}\) \(F = \frac{MS_A}{MS_{AB}}\)
Factor B \(SS_B\) \(b-1\) \(MS_B = \frac{SS_B}{b-1}\) \(F = \frac{MS_B}{MS_{AB}}\)
Interacció AB \(SS_{AB}\) \((a-1)(b-1)\) \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) \(F = \frac{MS_{AB}}{MS_E}\)
Error \(SS_E\) \(ab(n-1)\) \(MS_E = \frac{SS_E}{ab(n-1)}\)
Total \(SS_T\) \(N-1\)

1.2.4 Resolució de l’exemple 3

res.aov<-aov(compt~prep*analist,data=dades)
anova(res.aov)
Df Sum Sq Mean Sq F value Pr(>F)
prep 2 3362.0000000 1681.0000000 395.5294118 0.0000004
analist 1 44.0833333 44.0833333 10.3725490 0.0181225
prep:analist 2 0.6666667 0.3333333 0.0784314 0.9254977
Residuals 6 25.5000000 4.2500000 NA NA

Els resultats de l’ANOVA encara no són correctes. R assumeix per defecte que tots els factors són fixos.

1.2.4.1 Modificació manual de la taula ANOVA

## Modificació de la taula a ma
taula<-anova(res.aov)
taula[1,4]<-taula[1,3]/taula[3,3]
taula[2,4]<-taula[2,3]/taula[3,3]
taula[1,5]<-1-pf(taula[1,4],taula[1,1],taula[3,1])
taula[2,5]<-1-pf(taula[2,4],taula[2,1],taula[3,1])
taula
Df Sum Sq Mean Sq F value Pr(>F)
prep 2 3362.0000000 1681.0000000 5.04300e+03 0.0001983
analist 1 44.0833333 44.0833333 1.32250e+02 0.0074767
prep:analist 2 0.6666667 0.3333333 7.84314e-02 0.9254977
Residuals 6 25.5000000 4.2500000 NA NA

1.2.5 Estimació dels paràmetres del model. Components de la variància.

  • \(\hat{\mu} = \overline{Y}_{...}\)

  • \(\hat{\sigma}^2 = MS_E\)

  • \(\hat{\sigma}^2_A = \frac{MS_A-MS_{AB}}{bn}\)

  • \(\hat{\sigma}^2_B = \frac{MS_B-MS_{AB}}{an}\)

  • \(\hat{\sigma}^2_{AB} = \frac{MS_{AB}-MS_E}{n}\)

Aquestes estimacions s’obtenen pel mètode dels moments a partir de les esperances dels quadrats mitjans

\[ E(MS_E) = \sigma^2 \,\,\,\, ; \,\,\,\, E(MS_A) = \sigma^2 + bn \sigma^2_A+ n \sigma^2_{AB} \] \[ E(MS_B) = \sigma^2 + an \sigma^2_B+ n \sigma^2_{AB} \,\,\,\, ; \,\,\,\, E(MS_{AB}) = \sigma^2 + n \sigma^2_{AB} \]

1.2.6 Estimació de components de la variància a l’exemple 3

taula
Df Sum Sq Mean Sq F value Pr(>F)
prep 2 3362.0000000 1681.0000000 5.04300e+03 0.0001983
analist 1 44.0833333 44.0833333 1.32250e+02 0.0074767
prep:analist 2 0.6666667 0.3333333 7.84314e-02 0.9254977
Residuals 6 25.5000000 4.2500000 NA NA
  • \(\hat{\sigma}^2 = MS_E =\) 4.25.

  • \(\hat{\sigma}^2_{AB} = \frac{MS_{AB}-MS_E}{n} =\) -1.9583333 \(\Rightarrow 0\).

  • \(\hat{\sigma}^2_A = \frac{MS_A-MS_{AB}}{bn} =\) 420.1666667.

  • \(\hat{\sigma}^2_B = \frac{MS_B-MS_{AB}}{an} =\) 7.2916667.

compA<-(taula[1,3]-taula[3,3])/(4)
compB<-(taula[2,3]-taula[3,3])/(6)
compAB<-(taula[3,3]-taula[4,3])/(2)
compE<-taula[4,3]
if (compA<0){compA<-0}
if (compAB<0){compAB<-0}
if (compB<0){compB<-0}
perA<-100*compA/(compA+compB+compAB+compE)
perB<-100*compB/(compA+compB+compAB+compE)
perAB<-100*compAB/(compA+compB+compAB+compE)
perE<-100*compE/(compA+compB+compAB+compE)
c(compA,compB,compAB,compE)
## [1] 420.166667   7.291667   0.000000   4.250000
c(perA,perB,perAB,perE)
## [1] 97.326513  1.689026  0.000000  0.984461

1.2.6.1 Resolució amb el package mixlm

library(mixlm)
Anova(lm(compt~r(prep)*r(analist),data=dades),type='III')
## Analysis of variance (unrestricted model)
## Response: compt
##              Mean Sq  Sum Sq Df F value Pr(>F)
## prep         1681.00 3362.00  2 5043.00 0.0002
## analist        44.08   44.08  1  132.25 0.0075
## prep:analist    0.33    0.67  2    0.08 0.9255
## Residuals       4.25   25.50  6       -      -
## 
##                Err.term(s) Err.df VC(SS)
## 1 prep                 (3)      2 420.17
## 2 analist              (3)      2   7.29
## 3 prep:analist         (4)      6  -1.96
## 4 Residuals              -      -   4.25
## (VC = variance component)
## 
##              Expected mean squares
## prep           (4) + 2 (3) + 4 (1)
## analist        (4) + 2 (3) + 6 (2)
## prep:analist   (4) + 2 (3)        
## Residuals      (4)

1.2.7 Resolució alternativa de l’exemple 3 amb el package lme4

Ho farem només amb el package lme4, nlme amb la funció lme() sempre assumeix que els factors estan jerarquitzats (niats) i és molt difícil modelar factors aleatoris creuats.

1.2.7.1 Package lme4

library(lme4)
library(lmerTest)
res.lmer<-lmer(compt~1+(1|prep)+(1|analist)+(1|prep:analist))
summary(res.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: compt ~ 1 + (1 | prep) + (1 | analist) + (1 | prep:analist)
## 
## REML criterion at convergence: 61.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3909 -0.6030  0.1357  0.5627  1.3738 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  prep:analist (Intercept)   0.000   0.000  
##  prep         (Intercept) 419.435  20.480  
##  analist      (Intercept)   6.802   2.608  
##  Residual                   3.271   1.809  
## Number of obs: 12, groups:  prep:analist, 6; prep, 3; analist, 2
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  341.750     11.979   2.095   28.53 0.000949 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
confint(res.lmer)
##                   2.5 %     97.5 %
## .sig01        0.0000000   3.426253
## .sig02        8.8870627  49.709260
## .sig03        0.6087447  18.734220
## .sigma        1.1866227   3.247714
## (Intercept) 314.4782199 369.021888
ranova(res.lmer)
npar logLik AIC LRT Df Pr(>Chisq)
5 -30.91114 71.82228 NA NA NA
(1 | prep) 4 -37.50172 83.00345 13.181166 1 0.0002828
(1 | analist) 4 -32.67240 73.34479 3.522511 1 0.0605407
(1 | prep:analist) 4 -30.91114 69.82228 0.000000 1 1.0000000

1.2.8 Dos factors aleatoris amb una observació per casella

En aquest cas no podem incloure la interacció en el model i el model que resta és

\[ Y_{ij} = \mu + A_i + B_j +\epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b \] on

  • \(\mu\) és la mitjana general,

  • \(A_i \sim N(0,\sigma_A)\),

  • \(B_j \sim N(0,\sigma_B)\),

  • \(\epsilon_{ij}\) és l’error aleatori corresponent a l’observació \(i,j\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ij} \sim N(0,\sigma)\).

  • Assumirem també que les v.a. \(A_i\), \(B_j\) i \(\epsilon_{ijk}\) són v.a. independents.

En aquest cas \(SS_T = SS_A + SS_B + SS_E\) adaptant les fórmules al cas \(n=1\).

Els únics contrasts són:

  • Variabilitat aportada pel factor A \[ H_0: \sigma_A^2 = 0 \] \[ H_1: \sigma_A^2 > 0\]

  • Variabilitat aportada pel factor B \[ H_0: \sigma_B^2 = 0 \] \[ H_1: \sigma_B^2 > 0\]

La taula ANOVA queda

Font de variació Suma de quadrats Graus de llibertat Quadrats mitjans Estadístic F
Factor A \(SS_A\) \(a-1\) \(MS_A = \frac{SS_A}{a-1}\) \(F = \frac{MS_A}{MS_E}\)
Factor B \(SS_B\) \(b-1\) \(MS_B = \frac{SS_B}{b-1}\) \(F = \frac{MS_B}{MS_E}\)
Error \(SS_E\) \((a-1)(b-1)\) \(MS_E = \frac{SS_E}{(a-1)(b-1)}\)
Total \(SS_T\) \(N-1\)

1.2.8.1 Estimacions dels paràmetres del model

  • \(\hat{\mu} = \overline{Y}_{...}\)

  • \(\hat{\sigma}^2 = MS_E\)

  • \(\hat{\sigma}^2_A = \frac{MS_A-MS_E}{b}\)

  • \(\hat{\sigma}^2_B = \frac{MS_B-MS_E}{a}\)

1.3 Un factor fix i un factor aleatori

1.3.1 Exemple 4. Estudi d’adhesió cel·lular

En un estudi d’adhesió cel·lular es van analitzar 4 tipus cel·lulars segons la seva capacitat d’adhesió a un determinat substrat. Per raons logístiques es van realitzar tres experiments independents en tres dies diferents. En cada experiment es van analitzar els 4 tipus cel·lulars amb un total de 5 rèpliques per cada tipus. Volem valorar la influència en la variable resposta dels diferents factors que intervenen a l’estudi.

Els resultats són els següents

adhesio1<-c(0.89,0.95,1.31,1.01,1.03,1.24,0.99,0.65,0.76,0.81,
            0.94,0.77,0.78,0.72,0.59,1.36,1.1,0.61,1.14,0.91)
adhesio2<-c(0.69,0.86,0.69,0.75,0.77,0.46,0.97,0.83,0.57,1.05,
            0.69,0.54,0.78,0.78,0.76,1.04,1.01,0.89,0.84,0.82)
adhesio3<-c(0.36,0.46,0.25,0.93,0.74,0.18,0.32,0.31,0.56,0.38,
            0.25,0.19,0.19,0.47,0.15,0.8,0.19,0.62,0.52,0.54)
adhesio<-c(adhesio1,adhesio2,adhesio3)
cell<-rep(gl(4,5),3)
experiment<- gl(3,20)
myData <- data.frame(experiment,cell,adhesio)

El factor important és el tipus cel·lular, factor fix amb 4 nivells, però també volem estudiar la variabilitat que pot introduir el dia en que es fa l’experiment, per tant incluïm el factor experiment en el model. Tenim per tant un factor fix i un factor aleatori.

interaction.plot(myData$experiment,myData$cell,response=myData$adhesio,main='Interaction Plot',col='red',trace.label = 'cell',xlab='Experiment',ylab='Adhesió')

1.3.2 Dissenys mixtos. Un factor fix i un factor aleatori

Suposem ara que tenim un experiment amb dues fonts de variabilitat controlades \(A\) i \(B\) però on una es fixa \(A\) i l’altre aleatòria \(B\), es a dir els seus nivells són una mostra aleatòria de mida \(b\) d’un conjunt més gran de possibles nivells.

El model lineal és ara

\[ Y_{ijk} = \mu + \alpha_i + B_j + (\alpha B)_{ij} +\epsilon_{ijk} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,n_{ij} \]

on

  • \(Y_{ijk}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\),

  • \(\mu\) és la mitjana general,

  • \(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor A,

  • \(B_j\) és una variable aleatòria, no un paràmetre constant. Assumim que \(B_j \sim N(0,\sigma_B)\),

  • \((\alpha B)_{ij}\) és una variable aleatòria, no un paràmetre constant. Assumim que \((\alpha B)_{ij} \sim N(0,\sigma_{AB})\). (Tota interacció on al menys un dels termes sigui aleatori és aleatòria).

Assumim també la restricció \(\sum_{i=1}^a \alpha_i = 0\).

  • \(\epsilon_{ijk}\) és l’error aleatori corresponent a l’observació \(i,j,k\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ijk} \sim N(0,\sigma)\).

  • Assumirem també que les v.a. \(B_j\), \((\alpha B)_{ij}\) i \(\epsilon_{ijk}\) són v.a. independents.

i per tant cada observació experimental

\[ Y_{ijk} \sim N(\mu+\alpha_i,\sqrt{\sigma^2_B+\sigma^2_{AB}+\sigma^2} )\]

1.3.2.1 Dependència entre les observacions

\[ Cov(Y_{ijk},Y_{i'j'k'}) = \left\{ \begin{array}{ll} \sigma^2_B+\sigma^2_{AB} & \mbox{ si } i=i', j=j', k \ne k' \\ \sigma^2_B & \mbox{ si } i \ne i', j = j' \\ 0 & \mbox{ si } i \ne i', j \ne j' \end{array} \right. \]

Evidentment si \(i=i', j=j', k = k'\) llavors

\[ Cov(Y_{ijk},Y_{ijk})=Var(Y_{ijk})=\sigma^2_B+\sigma^2_{AB}+\sigma^2\]

1.3.3 Dissenys mixtos. Contrasts d’hipòtesis

Els contrasts de més interès ara són els següents

  • Hipòtesi sobre la significació del factor A \[ H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0 \] \[ H_1: \alpha_i \ne \alpha_{i'} \mbox{ per algun } i \ne i' \,\,\, ( \mbox{ equivalentment }H_1: \exists \alpha_i \ne 0)\]

  • Hipòtesi sobre la variabilitat aportada pel factor B \[ H_0: \sigma_B^2 = 0 \] \[ H_1: \sigma_B^2 > 0\]

  • Hipòtesi sobre la variabilitat aportada per la interacció \[ H_0: \sigma_{AB}^2 = 0 \] \[ H_1: \sigma_{AB}^2 > 0\]

Les esperances dels quadrats mitjans són les següents

  • \(E(MS_E) = \sigma^2\)
  • \(E(MS_A) = \sigma^2 + n \sigma^2_{AB}+ \frac{ bn \sum_{i=1}^a \alpha_i^2}{a-1}\)
  • \(E(MS_B) = \sigma^2 + an \sigma^2_B+ n \sigma^2_{AB}\)
  • \(E(MS_{AB}) = \sigma^2 + n \sigma^2_{AB}\)

Per tant els estadístics F per decidir en cadascú dels contrasts anteriors són

\[ F = \frac{MS_A}{MS_{AB}} \sim F(a-1,(a-1)(b-1)) \]

\[ F = \frac{MS_B}{MS_{AB}} \sim F(b-1,(a-1)(b-1)) \]

\[ F = \frac{MS_{AB}}{MS_E} \sim F((a-1)(b-1),ab(n-1)) \]

i la taula ANOVA és la següent:

Font de variació Suma de quadrats Graus de llibertat Quadrats mitjans Estadístic F
Factor A \(SS_A\) \(a-1\) \(MS_A = \frac{SS_A}{a-1}\) \(F = \frac{MS_A}{MS_{AB}}\)
Factor B \(SS_B\) \(b-1\) \(MS_B = \frac{SS_B}{b-1}\) \(F = \frac{MS_B}{MS_{AB}}\)
Interacció AB \(SS_{AB}\) \((a-1)(b-1)\) \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) \(F = \frac{MS_{AB}}{MS_E}\)
Error \(SS_E\) \(ab(n-1)\) \(MS_E = \frac{SS_E}{ab(n-1)}\)
Total \(SS_T\) \(N-1\)

1.3.3.1 Estimació dels paràmetres del model. Components de la variància.

  • \(\hat{\mu} = \overline{Y}_{...}\)

  • \(\hat{\alpha_i} = \overline{Y}_{i..} - \overline{Y}_{...} \,\,\,\, (\alpha_i = \mu_i - \mu)\)

  • \(\hat{\sigma}^2 = MS_E\)

  • \(\hat{\sigma}^2_B = \frac{MS_B-MS_{AB}}{an}\)

  • \(\hat{\sigma}^2_{AB} = \frac{MS_{AB}-MS_E}{n}\)

1.3.4 Dissenys mixtos. Model restringit

Suposem ara que tenim un experiment amb dues fonts de variabilitat controlades \(A\) i \(B\) però on una es fixa \(A\) i l’altre aleatòria \(B\), es a dir els seus nivells són una mostra aleatòria de mida \(b\) d’un conjunt més gran de possibles nivells.

El model lineal és

\[ Y_{ijk} = \mu + \alpha_i + B_j + (\alpha B)_{ij} +\epsilon_{ijk} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,n_{ij} \]

Però ara modifiquem les restriccions als paràmetres on

  • \(Y_{ijk}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\),

  • \(\mu\) és una mitjana general,

  • \(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor A amb \(\sum_{i=1}^a \alpha_i = 0\),

  • \(B_j \sim N(0,\sigma_B)\),

  • \(\epsilon_{ijk}\) és l’error aleatori corresponent a l’observació \(i,j,k\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ijk} \sim N(0,\sigma)\).

La nova restricció s’imposa sobre \((\alpha B)_{ij}\):

\[ (\alpha B)_{ij} \sim N \left( 0 , \sqrt{\frac{a-1}{a}}\sigma_{AB} \right) \,\,\,\, , \,\,\,\, \sum_{i=1}^a (\alpha B)_{ij} = 0, \,\,\, j=1,...,b \]

A causa d’aquesta nova restricció el model es coneix com model restringit. La forma de la variància de \((\alpha B)_{ij}\) simplifica algunes expressions de les esperances de les sumes de quadrats.

1.3.5 Model restringit. Contrasts d’hipòtesis

  • Hipòtesi sobre la significació del factor A \[ H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0 \] \[ H_1: \alpha_i \ne \alpha_{i'} \mbox{ per algun } i \ne i' \,\,\, ( \mbox{ equivalentment }H_1: \exists \alpha_i \ne 0)\]

  • Hipòtesi sobre la variabilitat aportada pel factor B \[ H_0: \sigma_B^2 = 0 \] \[ H_1: \sigma_B^2 > 0\]

  • Hipòtesi sobre la variabilitat aportada per la interacció \[ H_0: \sigma_{AB}^2 = 0 \] \[ H_1: \sigma_{AB}^2 > 0\]

Les esperances dels quadrats mitjans són les següents

  • \(E(MS_E) = \sigma^2\)
  • \(E(MS_A) = \sigma^2 + n \sigma^2_{AB}+ \frac{ bn \sum_{i=1}^a \alpha_i^2}{a-1}\)
  • \(E(MS_B) = \sigma^2 + an \sigma^2_B\)
  • \(E(MS_{AB}) = \sigma^2 + n \sigma^2_{AB}\)

Per tant els estadístics F per decidir en cadascú dels contrasts anteriors són

\[ F = \frac{MS_A}{MS_{AB}} \sim F(a-1,(a-1)(b-1)) \]

\[ F = \frac{MS_B}{MS_{E}} \sim F(b-1,ab(n-1)) \]

\[ F = \frac{MS_{AB}}{MS_E} \sim F((a-1)(b-1),ab(n-1)) \]

i la taula ANOVA és la següent

Font de variació Suma de quadrats Graus de llibertat Quadrats mitjans Estadístic F
Factor A \(SS_A\) \(a-1\) \(MS_A = \frac{SS_A}{a-1}\) \(F = \frac{MS_A}{MS_{AB}}\)
Factor B \(SS_B\) \(b-1\) \(MS_B = \frac{SS_B}{b-1}\) \(F = \frac{MS_B}{MS_{E}}\)
Interacció AB \(SS_{AB}\) \((a-1)(b-1)\) \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) \(F = \frac{MS_{AB}}{MS_E}\)
Error \(SS_E\) \(ab(n-1)\) \(MS_E = \frac{SS_E}{ab(n-1)}\)
Total \(SS_T\) \(N-1\)

1.3.5.1 Estimació dels paràmetres del model. Components de la variància.

  • \(\hat{\mu} = \overline{Y}_{...}\)

  • \(\hat{\alpha_i} = \overline{Y}_{i..} - \overline{Y}_{...} \,\,\,\, (\alpha_i = \mu_i - \mu)\)

  • \(\hat{\sigma}^2 = MS_E\)

  • \(\hat{\sigma}^2_B = \frac{MS_B-MS_{E}}{an}\)

  • \(\hat{\sigma}^2_{AB} = \frac{MS_{AB}-MS_E}{n}\)

1.3.6 Restringit vs. no restringit. Quin model utilitzar?

  • Algunes persones i alguns paquets de software tenen preferències per una opció per defecte i l’apliquen sempre.

  • Però ens podem preguntar si hi han situacions on un model s’ajusta millor que l’altre.

  • La conseqüència més immediata de la restricció imposada al model restringit és que

\[ Cov ((\alpha B)_{ij},(\alpha B)_{i'j}) = -\sigma^2_{AB}/a. \]

Es a dir els termes d’interacció dintre de cada nivell del factor aleatori estan correlacionats negativament.

Això es tradueix en la següent estructura de dependència entre observacions pertanyents al mateix nivell del factor aleatori en el model restringit:

\[ Cov (Y_{ijk},Y_{i'jk'}) = \sigma_B^2-\frac{\sigma_{AB}^2}{a}. \]

Mentre que en el model no restringit

\[ Cov (Y_{ijk},Y_{i'jk'}) = \sigma_B^2. \]

Com a detall important en el model restringit la covariància pot ser negativa o positiva mentre que en el model no restringit la covariància sempre és positiva.

  • Una situació en que això es reflecteix és en experiments amb recursos limitats.

Per exemple imaginem que estem comparant el creixement de dues espècies de plantes (factor fix) i que la parcel·la on creixen és un factor aleatori de b nivells al tractar-se de parcel·les seleccionades aleatòriament. Cada parcel·la es divideix en 2n unitats experimentals i es produeix una assignació aleatòria dels individus de forma que cada espècie té n individus en cada parcel·la. Degut a que els recursos (nutrients, aigua, espai) a cada parcel·la son limitats, esperem correlacions negatives entre els pesos finals de les plantes de cada parcel·la. Per tant en aquest cas semblaria millor el model restringit.

1.3.7 Resolució de l’exemple 4

mixed.aov<-aov(adhesio~cell*experiment,data=myData)
anova(mixed.aov)
Df Sum Sq Mean Sq F value Pr(>F)
cell 3 0.5753467 0.1917822 5.3827311 0.0028246
experiment 2 2.7526633 1.3763317 38.6293533 0.0000000
cell:experiment 6 0.1201633 0.0200272 0.5621019 0.7581959
Residuals 48 1.7102000 0.0356292 NA NA

1.3.7.1 Resolució manual del model no restringit

taula<-anova(mixed.aov)
taula[1,4]<-taula[1,3]/taula[3,3]
taula[2,4]<-taula[2,3]/taula[3,3]
taula[1,5]<-1-pf(taula[1,4],taula[1,1],taula[3,1])
taula[2,5]<-1-pf(taula[2,4],taula[2,1],taula[3,1])
taula
Df Sum Sq Mean Sq F value Pr(>F)
cell 3 0.5753467 0.1917822 9.5760770 0.0105230
experiment 2 2.7526633 1.3763317 68.7230436 0.0000732
cell:experiment 6 0.1201633 0.0200272 0.5621019 0.7581959
Residuals 48 1.7102000 0.0356292 NA NA

1.3.7.2 Resolució amb el package mixlm

Anova(lm(adhesio~cell*r(experiment),data=myData,unrestricted=T),type='III')
## Analysis of variance (unrestricted model)
## Response: adhesio
##                 Mean Sq Sum Sq Df F value Pr(>F)
## cell               0.19   0.58  3    9.58 0.0105
## experiment         1.38   2.75  2   68.72 0.0001
## cell:experiment    0.02   0.12  6    0.56 0.7582
## Residuals          0.04   1.71 48       -      -
## 
##                   Err.term(s) Err.df   VC(SS)
## 1 cell                    (3)      6    fixed
## 2 experiment              (3)      6  0.06782
## 3 cell:experiment         (4)     48 -0.00312
## 4 Residuals                 -      -  0.03563
## (VC = variance component)
## 
##                 Expected mean squares
## cell            (4) + 5 (3) + 15 Q[1]
## experiment      (4) + 5 (3) + 20 (2) 
## cell:experiment (4) + 5 (3)          
## Residuals       (4)

1.3.7.3 Resolució manual del model restringit

taula<-anova(mixed.aov)
taula[1,4]<-taula[1,3]/taula[3,3]
taula[1,5]<-1-pf(taula[1,4],taula[1,1],taula[3,1])
taula
Df Sum Sq Mean Sq F value Pr(>F)
cell 3 0.5753467 0.1917822 9.5760770 0.0105230
experiment 2 2.7526633 1.3763317 38.6293533 0.0000000
cell:experiment 6 0.1201633 0.0200272 0.5621019 0.7581959
Residuals 48 1.7102000 0.0356292 NA NA

1.3.7.4 Resolució amb el package mixlm

Anova(lm(adhesio~cell*r(experiment),data=myData,unrestricted=F),type='III')
## Analysis of variance (restricted model)
## Response: adhesio
##                 Mean Sq Sum Sq Df F value Pr(>F)
## cell               0.19   0.58  3    9.58 0.0105
## experiment         1.38   2.75  2   38.63 0.0000
## cell:experiment    0.02   0.12  6    0.56 0.7582
## Residuals          0.04   1.71 48       -      -
## 
##                   Err.term(s) Err.df   VC(SS)
## 1 cell                    (3)      6    fixed
## 2 experiment              (4)     48  0.06704
## 3 cell:experiment         (4)     48 -0.00312
## 4 Residuals                 -      -  0.03563
## (VC = variance component)
## 
##                 Expected mean squares
## cell            (4) + 5 (3) + 15 Q[1]
## experiment      (4) + 20 (2)         
## cell:experiment (4) + 5 (3)          
## Residuals       (4)

1.3.8 Comparacions múltiples

L’única part complicada dels contrasts o comparacions per parelles sobre els factors principals fixos és assegurar-se que s’utilitza el terme d’error correcte si el model conté factors aleatoris. En aquestes situacions, el terme d’error per al factor fix, i per tant per qualsevol contrast sobre el valor mitjà d’aquests factors, sol ser la interacció en lloc del residual.

df<-taula[3,1]
MSerror<-taula[3,3]
HSD.test(myData$adhesio,myData$cell,df,MSerror,console=T,group=F)
## 
## Study: myData$adhesio ~ myData$cell
## 
## HSD Test for myData$adhesio 
## 
## Mean Square Error:  0.02002722 
## 
## myData$cell,  means
## 
##   myData.adhesio       std  r         se  Min  Max   Q25  Q50   Q75
## 1      0.7793333 0.2731684 15 0.03653968 0.25 1.31 0.690 0.77 0.940
## 2      0.6720000 0.3114299 15 0.03653968 0.18 1.24 0.420 0.65 0.900
## 3      0.5733333 0.2619887 15 0.03653968 0.15 0.94 0.360 0.69 0.775
## 4      0.8260000 0.2951465 15 0.03653968 0.19 1.36 0.615 0.84 1.025
## 
## Alpha: 0.05 ; DF Error: 6 
## Critical Value of Studentized Range: 4.895599 
## 
## Comparison between treatments means
## 
##        difference pvalue signif.         LCL         UCL
## 1 - 2  0.10733333 0.2605         -0.07155029  0.28621696
## 1 - 3  0.20600000 0.0277       *  0.02711638  0.38488362
## 1 - 4 -0.04666667 0.8043         -0.22555029  0.13221696
## 2 - 3  0.09866667 0.3163         -0.08021696  0.27755029
## 2 - 4 -0.15400000 0.0881       . -0.33288362  0.02488362
## 3 - 4 -0.25266667 0.0109       * -0.43155029 -0.07378304

1.3.9 Un factor fix i un aleatori amb lmer(lme4) (Opcional)

Per defecte lmer() fa servir el model no restringit. No és possible definir lmer() amb el model restringit.

library(lme4)
library(lmerTest)
res.lmer<-lmer(adhesio~cell+(1|experiment)+(1|cell:experiment),data=myData)
summary(res.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: adhesio ~ cell + (1 | experiment) + (1 | cell:experiment)
##    Data: myData
## 
## REML criterion at convergence: -12.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3140 -0.5718 -0.0934  0.5960  2.3662 
## 
## Random effects:
##  Groups          Name        Variance Std.Dev.
##  cell:experiment (Intercept) 0.00000  0.0000  
##  experiment      (Intercept) 0.06712  0.2591  
##  Residual                    0.03390  0.1841  
## Number of obs: 60, groups:  cell:experiment, 12; experiment, 3
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)  0.77933    0.15695  2.30615   4.966   0.0284 * 
## cell2       -0.10733    0.06723 53.99995  -1.597   0.1162   
## cell3       -0.20600    0.06723 53.99995  -3.064   0.0034 **
## cell4        0.04667    0.06723 53.99995   0.694   0.4906   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) cell2  cell3 
## cell2 -0.214              
## cell3 -0.214  0.500       
## cell4 -0.214  0.500  0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
confint(res.lmer)
##                   2.5 %      97.5 %
## .sig01       0.00000000  0.07439780
## .sig02       0.10643502  0.63124876
## .sigma       0.15073129  0.21789385
## (Intercept)  0.42904003  1.12962474
## cell2       -0.23781143  0.02314497
## cell3       -0.33647809 -0.07552170
## cell4       -0.08381143  0.17714497
anova(res.lmer)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
cell 0.5753467 0.1917822 3 53.99995 5.658015 0.001911
ranova(res.lmer)
npar logLik AIC LRT Df Pr(>Chisq)
7 6.1845980 1.630804 NA NA NA
(1 | experiment) 6 -0.7499014 13.499803 13.869 1 0.000196
(1 | cell:experiment) 6 6.1845980 -0.369196 0.000 1 1.000000

1.3.9.1 Comparacions múltiples

Per les comparacions múltiples després d’un ajust amb lmer tenim diverses alternatives, presentarem la resolució amb el package emmeans

library(emmeans)
emm<-emmeans(res.lmer,'cell')
pairs(emm,adjust='tukey')
##  contrast      estimate     SE df t.ratio p.value
##  cell1 - cell2   0.1073 0.0672  6   1.597  0.4456
##  cell1 - cell3   0.2060 0.0672  6   3.064  0.0797
##  cell1 - cell4  -0.0467 0.0672  6  -0.694  0.8958
##  cell2 - cell3   0.0987 0.0672  6   1.468  0.5080
##  cell2 - cell4  -0.1540 0.0672  6  -2.291  0.2021
##  cell3 - cell4  -0.2527 0.0672  6  -3.758  0.0357
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates
pwpp(emm)

## En forma equivalent
## emmeans(res.lmer, pairwise ~ cell, adjust = c("tukey"))

1.3.10 Un factor fix i un aleatori amb una observació per casella

En aquest cas no podem incloure la interacció en el model i el model lineal que resta és

\[ Y_{ij} = \mu + \alpha_i + B_j + \epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b \]

on

  • \(Y_{ij}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\),

  • \(\mu\) és la mitjana general,

  • \(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor A,

  • \(B_j\) és una variable aleatòria, no un paràmetre constant. Assumim que \(B_j \sim N(0,\sigma_B)\).

Assumim també la restricció \(\sum_{i=1}^a \alpha_i = 0\).

  • \(\epsilon_{ij}\) és l’error aleatori corresponent a l’observació \(i,j\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ij} \sim N(0,\sigma)\).

  • Assumirem també que les v.a. \(B_j\), \((\alpha B)_{ij}\) i \(\epsilon_{ij}\) són v.a. independents.

Font de variació Suma de quadrats Graus de llibertat Quadrats mitjans Estadístic F
Factor A \(SS_A\) \(a-1\) \(MS_A = \frac{SS_A}{a-1}\) \(F = \frac{MS_A}{MS_E}\)
Factor B \(SS_B\) \(b-1\) \(MS_B = \frac{SS_B}{b-1}\) \(F = \frac{MS_B}{MS_E}\)
Error \(SS_E\) \((a-1)(b-1)\) \(MS_E = \frac{SS_E}{(a-1)(b-1)}\)
Total \(SS_T\) \(N-1\)

1.3.10.1 Estimació dels paràmetres del model. Components de la variància.

  • \(\hat{\mu} = \overline{Y}_{..}\)

  • \(\hat{\alpha_i} = \overline{Y}_{i.} - \overline{Y}_{..} \,\,\,\, (\alpha_i = \mu_i - \mu)\)

  • \(\hat{\sigma}^2 = MS_E\)

  • \(\hat{\sigma}^2_B = \frac{MS_B-MS_{E}}{a}\)

1.3.10.2 Resolució amb mixlm (Opcional)

El paràmetre unrestricted no és necessari donat que en aquest model no existeix versió restringida.

myData1<-myData[c(1,6,11,16,21,26,31,36,41,46,51,56),]
library(mixlm)
myData1.lm<-lm(adhesio~cell+r(experiment),unrestricted=T,data=myData1)
Anova(myData1.lm,type='III')
## Analysis of variance (unrestricted model)
## Response: adhesio
##            Mean Sq Sum Sq Df F value Pr(>F)
## cell          0.14   0.42  3    6.15 0.0292
## experiment    0.51   1.01  2   22.03 0.0017
## Residuals     0.02   0.14  6       -      -
## 
##              Err.term(s) Err.df VC(SS)
## 1 cell               (3)      6  fixed
## 2 experiment         (3)      6 0.1206
## 3 Residuals            -      - 0.0229
## (VC = variance component)
## 
##            Expected mean squares
## cell                (3) + 3 Q[1]
## experiment          (3) + 4 (2) 
## Residuals           (3)

1.4 Dissenys amb blocs

  • A tot experiment pot haver un factor de soroll que pugui afectar la variabilitat dels resultats però en el que no estem directament interessats.

  • A vegades aquest factor és desconegut i no controlable, llavors l’aleatorització és l’eina que ens permet mitigar l’efecte d’aquest factor que queda incorporat a l’error experimental. De fet les pròpies unitats experimentals són una font de variabilitat.

  • Altres vegades és un factor conegut i controlable llavors podem incorporar-lo al disseny i eliminar de l’error experimental la variabilitat deguda a ell. Direm que tenim un factor bloc. No estem interessats de forma específica en contrastar el factor bloc, la raó d’incorporar-lo al disseny és eliminar la variabilitat atribuïda a ell de la variabilitat residual i incrementar la potència del test sobre el factor principal.

  • Construir un disseny tenint en compte possibles efectes bloc rep el nom en anglés de blocking i és un dels principis formulats per Ronald Fisher juntament amb l’aleatorització i la replicació.

“Block what you can; randomize what you cannot.”

  • Habitualment no són de interès les interaccions en que participa un factor bloc.

1.4.1 Disseny en blocs aleatoritzats (RCBD)

  • Un dels dissenys més utilitzats en aquest cas és el disseny en blocs aleatoritzats: S’anomena disseny en blocs aleatoritzats un experiment en el qual es pren el mateix nombre de mesures per tractament dins de cada bloc i l’ordre de les mesures dins de cada bloc es decideix aleatòriament.

1.4.2 Exemple 5: Influència de l’exercici físic

Es considera que amb l’exercici físic el pH de la sang tendeix a baixar, a causa de l’acumulació de \(CO_2\) i d’ions \(H^+\) a la musculatura. Que el pH sanguini baixi per sota de 6,8 (o que pugi per sobre de 7,4) pot provocar una crisi greu, amb possibilitat de mort. Sortosament hi ha mecanismes de regulació que tendeixen a mantenir el pH en nivells adequats. En un estudi als Estats Units es volia estudiar la relació entre el pH mitjà de la sang i l’exercici físic. Es consideraren tres règims d’exercici intens A, B i C, i per qüestions administratives, els individus del estudi foren triats al atzar dintre de 6 hospitals triats al atzar. Els pH sanguinis mesurats a l’experiment foren els següents:

Règim exercici Hosp 1 Hosp 2 Hosp 3 Hosp 4 Hosp 5 Hosp 6
A 7,29 7,20 7,15 7,23 7,2 7,1
B 7,26 7,16 7,08 7,11 7,16 6,89
C 7,18 7,13 6,99 7,03 7,11 6,90
pH    <- c( 7.29,7.20,7.15,7.23,7.20,7.10,7.26,7.16,7.08,
            7.11,7.16,6.89,7.18,7.13,6.99,7.03,7.11,6.90)
exer  <- factor(c(rep(1,6),rep(2,6),rep(3,6)))
hosp  <- factor(c(rep(1:6,3)))
pHdata  <- data.frame(exer,hosp,pH)
boxplot(pH~exer,pHdata,col='orange')

1.4.3 RCBD: Model

1.4.3.1 Bloc Fix

\[ Y_{ij} = \mu + \alpha_i + \beta_j +\epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b \]

  • \(\mu\) és la mitjana general,

  • \(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor tractament,

  • \(\beta_j\) és un paràmetre del model que mesura l’efecte del nivell \(j\) del factor de soroll (bloc),

  • \(\epsilon_{ij}\) és l’error aleatori corresponent a l’observació \(i,j\) i dels quals suposem que són independents i idènticament distribuïts amb \(E(\epsilon_{ij}) = 0\) i \(Var(\epsilon_{ij})=\sigma^2\).

La resolució és idèntica al cas de dos factors amb una rèplica per casella vist anteriorment.

\(SS_T = SS_A + SS_B + SS_E\) adaptant les fórmules al cas \(n=1\).

Tinguem en compte que l’única hipòtesi d’interès és sobre el factor primari, el factor tractament, no ens interessa el contrast sobre el bloc tot i que es podria fer.

\[ H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0 \] \[ H_1: \alpha_i \ne \alpha_{i'} \mbox{ per algun } i \ne i' \,\,\, ( \mbox{ equivalentment }H_1: \exists \alpha_i \ne 0)\]

Podem utilitzar una versió simplificada de la taula ANOVA

Font de variació Suma de quadrats Graus de llibertat Quadrats mitjans Estadístic F
Factor A \(SS_A\) \(a-1\) \(MS_A = \frac{SS_A}{a-1}\) \(F = \frac{MS_A}{MS_E}\)
Factor B \(SS_B\) \(b-1\) \(MS_B = \frac{SS_B}{b-1}\)
Error \(SS_E\) \((a-1)(b-1)\) \(MS_E = \frac{SS_E}{(a-1)(b-1)}\)
Total \(SS_T\) \(N-1\)

1.4.3.2 Bloc aleatori

És molt habitual que tal i com passa a l’exemple 5 el factor bloc sigui un factor aleatori. En aquest cas la taula i les consideracions són idèntiques. Tal sols cal tenir en compte que si es vol contrastar l’efecte del bloc cal interpretar de forma adequada la hipòtesis nul·la corresponent i adaptar els paràmetres del model.

\[ Y_{ij} = \mu + \alpha_i + B_j +\epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b \]

  • \(\mu\) és la mitjana general,

  • \(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor tractament,

  • \(B_j\) és una variable aleatòria \(B_j \sim N(0,\sigma_B)\) que mesura la variabilitat aportada pel factor de soroll (bloc),

  • \(\epsilon_{ij}\) és l’error aleatori corresponent a l’observació \(i,j\) i dels quals suposem que són independents i idènticament distribuïts amb \(E(\epsilon_{ij}) = 0\) i \(Var(\epsilon_{ij})=\sigma^2\). Són independents de les \(B_j\).

  • La resolució és idèntica perquè no ens interessa contrastar el factor bloc.

1.4.4 Resolució de l’exemple 5

## Models amb blocs

pHdata.aov <- aov(pH ~ exer + hosp, pHdata)           
anova(pHdata.aov)
Df Sum Sq Mean Sq F value Pr(>F)
exer 2 0.0584111 0.0292056 18.49754 0.0004363
hosp 5 0.1354944 0.0270989 17.16327 0.0001267
Residuals 10 0.0157889 0.0015789 NA NA
## Comparacions multiples
library(agricolae)
HSD.test(pHdata.aov, "exer",console=T) 
## 
## Study: pHdata.aov ~ "exer"
## 
## HSD Test for pH 
## 
## Mean Square Error:  0.001578889 
## 
## exer,  means
## 
##         pH        std r         se  Min  Max    Q25   Q50    Q75
## 1 7.195000 0.06534524 6 0.01622184 7.10 7.29 7.1625 7.200 7.2225
## 2 7.110000 0.12393547 6 0.01622184 6.89 7.26 7.0875 7.135 7.1600
## 3 7.056667 0.10308572 6 0.01622184 6.90 7.18 7.0000 7.070 7.1250
## 
## Alpha: 0.05 ; DF Error: 10 
## Critical Value of Studentized Range: 3.876777 
## 
## Minimun Significant Difference: 0.06288846 
## 
## Treatments with the same letter are not significantly different.
## 
##         pH groups
## 1 7.195000      a
## 2 7.110000      b
## 3 7.056667      b

Anem a veure que succeeix si no incloem el factor hospital

## Model sense efecte bloc

pHdata2.aov<-aov(pH~exer,pHdata)
anova(pHdata2.aov)
Df Sum Sq Mean Sq F value Pr(>F)
exer 2 0.0584111 0.0292056 2.89578 0.0864035
Residuals 15 0.1512833 0.0100856 NA NA

1.4.4.1 Consideracions

  • El model assumeix que no hi ha interacció. Si la suposició és falsa és perillós de cara a interpretar els resultats. veure la consideració anterior sobre la prova d’aditivitat de Tukey.

  • Es un disseny molt eficient respecte al número de casos necessaris.

  • Són vàlids tots els mètodes de comparacions múltiples.

  • L’estimació de paràmetres és idèntica a la vista anteriorment.

  • Cal verificar les suposicions per poder aplicar l’ANOVA.

1.4.5 Resolució de l’exemple 5 amb els packages nlme i lme4

1.4.5.1 Package nlme

library(nlme)
pHdata.lme<-lme(pH~exer,data=pHdata,random=~1|hosp)
summary(pHdata.lme)
## Linear mixed-effects model fit by REML
##   Data: pHdata 
##         AIC       BIC   logLik
##   -24.60822 -21.06797 17.30411
## 
## Random effects:
##  Formula: ~1 | hosp
##         (Intercept)   Residual
## StdDev:  0.09223159 0.03973524
## 
## Fixed effects:  pH ~ exer 
##                 Value  Std.Error DF   t-value p-value
## (Intercept)  7.195000 0.04099910 10 175.49167  0.0000
## exer2       -0.085000 0.02294115 10  -3.70513  0.0041
## exer3       -0.138333 0.02294115 10  -6.02992  0.0001
##  Correlation: 
##       (Intr) exer2
## exer2 -0.28       
## exer3 -0.28   0.50
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.81043774 -0.54870477  0.09034733  0.46540220  1.33538484 
## 
## Number of Observations: 18
## Number of Groups: 6
intervals(pHdata.lme)
## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                  lower       est.       upper
## (Intercept)  7.1036483  7.1950000  7.28635168
## exer2       -0.1361161 -0.0850000 -0.03388394
## exer3       -0.1894494 -0.1383333 -0.08721727
## 
##  Random Effects:
##   Level: hosp 
##                      lower       est.     upper
## sd((Intercept)) 0.04773274 0.09223159 0.1782145
## 
##  Within-group standard error:
##      lower       est.      upper 
## 0.02563547 0.03973524 0.06159001
anova(pHdata.lme)
numDF denDF F-value p-value
(Intercept) 1 10 33678.19389 0.0000000
exer 2 10 18.49754 0.0004363

1.4.5.2 Package lme4

library(lme4)
pHdata.lmer<-lmer(pH~exer+(1|hosp),data=pHdata)
summary(pHdata.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pH ~ exer + (1 | hosp)
##    Data: pHdata
## 
## REML criterion at convergence: -34.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.81044 -0.54870  0.09035  0.46540  1.33538 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  hosp     (Intercept) 0.008507 0.09223 
##  Residual             0.001579 0.03974 
## Number of obs: 18, groups:  hosp, 6
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  7.19500    0.04100  6.19114 175.492 1.11e-12 ***
## exer2       -0.08500    0.02294 10.00000  -3.705 0.004074 ** 
## exer3       -0.13833    0.02294 10.00000  -6.030 0.000127 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) exer2 
## exer2 -0.280       
## exer3 -0.280  0.500
confint(pHdata.lmer)
##                   2.5 %      97.5 %
## .sig01       0.04932455  0.17271691
## .sigma       0.02548100  0.05753568
## (Intercept)  7.11029296  7.27970703
## exer2       -0.12956104 -0.04043896
## exer3       -0.18289437 -0.09377229
anova(pHdata.lmer)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
exer 0.0584111 0.0292056 2 10 18.49754 0.0004363
library(lmerTest)
ranova(pHdata.lmer)
npar logLik AIC LRT Df Pr(>Chisq)
5 17.304109 -24.60822 NA NA NA
(1 | hosp) 4 9.404553 -10.80911 15.79911 1 7.04e-05

1.4.6 Resolució de l’exemple 4 com un model on experiment és un bloc

Considerearem que el factor experiment de l’exemple 4 és un bloc i no inclourem la interacció amb el factor fix cell.

1.4.6.1 Resolució amb aov

ex4.bloc<-aov(adhesio~cell+experiment,data=myData)
anova(ex4.bloc)
Df Sum Sq Mean Sq F value Pr(>F)
cell 3 0.5753467 0.1917822 5.658024 0.001911
experiment 2 2.7526633 1.3763317 40.605004 0.000000
Residuals 54 1.8303633 0.0338956 NA NA
HSD.test(ex4.bloc,'cell',console=T,group=F)
## 
## Study: ex4.bloc ~ "cell"
## 
## HSD Test for adhesio 
## 
## Mean Square Error:  0.03389562 
## 
## cell,  means
## 
##     adhesio       std  r         se  Min  Max   Q25  Q50   Q75
## 1 0.7793333 0.2731684 15 0.04753638 0.25 1.31 0.690 0.77 0.940
## 2 0.6720000 0.3114299 15 0.04753638 0.18 1.24 0.420 0.65 0.900
## 3 0.5733333 0.2619887 15 0.04753638 0.15 0.94 0.360 0.69 0.775
## 4 0.8260000 0.2951465 15 0.04753638 0.19 1.36 0.615 0.84 1.025
## 
## Alpha: 0.05 ; DF Error: 54 
## Critical Value of Studentized Range: 3.748904 
## 
## Comparison between treatments means
## 
##        difference pvalue signif.         LCL         UCL
## 1 - 2  0.10733333 0.3891         -0.07087601  0.28554267
## 1 - 3  0.20600000 0.0174       *  0.02779066  0.38420934
## 1 - 4 -0.04666667 0.8989         -0.22487601  0.13154267
## 2 - 3  0.09866667 0.4638         -0.07954267  0.27687601
## 2 - 4 -0.15400000 0.1129         -0.33220934  0.02420934
## 3 - 4 -0.25266667 0.0023      ** -0.43087601 -0.07445733

1.4.6.2 Resolució amb el package mixlm

Anova(lm(adhesio~cell+r(experiment),data=myData),type='III')
## Analysis of variance (unrestricted model)
## Response: adhesio
##            Mean Sq Sum Sq Df F value Pr(>F)
## cell          0.19   0.58  3    5.66 0.0019
## experiment    1.38   2.75  2   40.61 0.0000
## Residuals     0.03   1.83 54       -      -
## 
##              Err.term(s) Err.df VC(SS)
## 1 cell               (3)     54  fixed
## 2 experiment         (3)     54 0.0671
## 3 Residuals            -      - 0.0339
## (VC = variance component)
## 
##            Expected mean squares
## cell               (3) + 15 Q[1]
## experiment         (3) + 20 (2) 
## Residuals          (3)

1.4.6.3 Resolució amb el package lme4

res.lmer<-lmer(adhesio~cell+(1|experiment),data=myData)
summary(res.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: adhesio ~ cell + (1 | experiment)
##    Data: myData
## 
## REML criterion at convergence: -12.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3140 -0.5718 -0.0934  0.5960  2.3662 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  experiment (Intercept) 0.06712  0.2591  
##  Residual               0.03390  0.1841  
## Number of obs: 60, groups:  experiment, 3
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)  0.77933    0.15695  2.30598   4.965   0.0284 * 
## cell2       -0.10733    0.06723 54.00000  -1.597   0.1162   
## cell3       -0.20600    0.06723 54.00000  -3.064   0.0034 **
## cell4        0.04667    0.06723 54.00000   0.694   0.4906   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr) cell2  cell3 
## cell2 -0.214              
## cell3 -0.214  0.500       
## cell4 -0.214  0.500  0.500
anova(res.lmer)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
cell 0.5753467 0.1917822 3 54 5.658024 0.001911
ranova(res.lmer)
npar logLik AIC LRT Df Pr(>Chisq)
6 6.184598 -0.369196 NA NA NA
(1 | experiment) 5 -16.179175 42.358350 44.72755 1 0
emm<-emmeans(res.lmer,'cell')
pairs(emm,adjust='tukey')
##  contrast      estimate     SE df t.ratio p.value
##  cell1 - cell2   0.1073 0.0672 54   1.597  0.3891
##  cell1 - cell3   0.2060 0.0672 54   3.064  0.0174
##  cell1 - cell4  -0.0467 0.0672 54  -0.694  0.8989
##  cell2 - cell3   0.0987 0.0672 54   1.468  0.4638
##  cell2 - cell4  -0.1540 0.0672 54  -2.291  0.1129
##  cell3 - cell4  -0.2527 0.0672 54  -3.758  0.0023
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 4 estimates
pwpp(emm)

1.5 El cas no balancejat

Eliminarem una observació del conjunt de dades del primer exemple i veurem què passa al fer un Anova intercanviant l’ordre dels factors.

immuno_1<-immuno[-1,]
anova(aov(anticos~vacuna*dosis,immuno_1))
Df Sum Sq Mean Sq F value Pr(>F)
vacuna 3 141.498359 47.1661197 11.3616797 0.0000236
dosis 2 8.178113 4.0890564 0.9849983 0.3835536
vacuna:dosis 6 5.375372 0.8958953 0.2158091 0.9692567
Residuals 35 145.296667 4.1513333 NA NA
anova(aov(anticos~dosis*vacuna,immuno_1))
Df Sum Sq Mean Sq F value Pr(>F)
dosis 2 7.902302 3.9511512 0.9517788 0.3958190
vacuna 3 141.774170 47.2580565 11.3838261 0.0000232
dosis:vacuna 6 5.375372 0.8958953 0.2158091 0.9692567
Residuals 35 145.296667 4.1513333 NA NA

Per què obtenim diferències en les taules Anova?

1.5.1 Dissenys balancejats vs no balancejats

Hi han moltes raons per preferir un disseny balancejat:

  • Els contrasts són molt més robusts a la suposició de normalitat i homoscedasticitat quan les mides mostrals són iguals.

  • L’estimació dels components de la variància és molt més complicada amb mides mostrals diferents.

Però més enllà d’aquestes raons en un disseny no balancejat no hi ha una única partició additiva de la suma de quadrats total. Això implica que hi han diferents formes de calcular les sumes de quadrats. En general s’anomenen Tipus I, II i III.

R per defecte utilitza el tipus I, a diferència d’altres programes com per exemple SPSS que utilitza el tipus III.

Fem servir la notació SS(A,B,AB) per indicar un model amb els factors A, B i la seva interacció AB; SS(A,B) el model sense la interacció; SS(B,AB) el model sense els efectes del factor A i així successivament. Llavors els diferents tipus de sumes de quadrats assignen els següents valors com a sumes de quadrats:

  • Tipus I (seqüencial) model A*B

    • SS(A) pel factor A
    • SS(B|A) = SS(A,B)-SS(A) pel factor B
    • SS(AB|B,A) = SS(A,B,AB)-SS(A,B) per a la interacció AB

    Els contrasts que es fan són: sobre el factor A, sobre el factor B després d’haver controlat pel factor A i sobre la interacció després d’haver controlat pels dos efectes principals.

  • Tipus II

    • SS(A|B) = SS(A,B)-SS(B) per el factor A
    • SS(B|A) = SS(A,B)-SS(A) per el factor B
  • Tipus III

    • SS(A|B,AB) = SS(A,B,AB)-SS(B,AB) per el factor A
    • SS(B|A,AB) = SS(A,B,AB)-SS(A,AB) per el factor B
    • SS(AB|B,A) = SS(A,B,AB)-SS(A,B) per a la interacció AB

La suma de quadrats de tipus I utilitza un procediment jeràrquic pel càlcul de les sumes de quadrats i l’ordre dels termes en el model té importància.

Molts autors recomanen utilitzar el tipus III que no està influenciat per l’ordre dels termes.

1.5.2 Comprovació de les fòrmules de les sumes de quadrats

Calculem manualment les sumes de quadrats dels diferents models

mod0<-lm(anticos~1,immuno_1)
SST<-sum(mod0$residuals^2)
modA<-lm(anticos~vacuna,immuno_1)
modB<-lm(anticos~dosis,immuno_1)
SSRA<-sum(modA$residuals^2)
SSA<-SST-SSRA
SSRB<-sum(modB$residuals^2)
SSB<-SST-SSRB
modA_B<-lm(anticos~vacuna+dosis,immuno_1)
SSRA_B<-sum(modA_B$residuals^2)
SSA_B<-SST-SSRA_B
modAB<-lm(anticos~vacuna*dosis,immuno_1)
SSRAB<-sum(modAB$residuals^2)
SSAB<-SST-SSRAB

X<-model.matrix(lm(anticos~vacuna*dosis,immuno_1))
interi<-X[,7:12]
modnB<-lm(anticos~vacuna+interi,immuno_1)
modnA<-lm(anticos~dosis+interi,immuno_1)
SSRnA<-sum(modnA$residuals^2)
SSRnB<-sum(modnB$residuals^2)
SSnA<-SST-SSRnA
SSnB<-SST-SSRnB
c(SST,SSA,SSB,SSA_B,SSAB,SSnA,SSnB)
## [1] 300.348511 141.498359   7.902302 149.676472 155.051844  12.099152 146.768099

Segons les taules anteriors les sumes de quadrats han de ser:

1.5.2.1 Tipus I amb vacuna*dosi

  • SS(vacuna): SSA = 141.4983591
  • SS(dosi): SSA_B - SSA = 8.1781128
  • SS(interacció): SSAB-SSA_B = 5.375372

1.5.2.2 Tipus I amb dosi*vacuna

  • SS(dosi): SSB = 7.9023023
  • SS(vacuna): SSA_B - SSB = 141.7741696
  • SS(interacció): SSAB - SSA_B =5.375372

1.5.2.3 Tipus II sense interacció

  • SS(vacuna): SSA_B - SSB = 141.7741696
  • SS(dosi): SSA_B - SSA = 8.1781128

1.5.2.4 Tipus III

  • SS(vacuna): SSAB - SSnA = 142.9526923
  • SS(dosi): SSAB - SSnB = 8.2837445
  • SS(interacció): SSAB - SSA-B = 5.375372

1.5.3 Suma de quadrats tipus III amb R

1.5.3.1 Amb el package sasLM

## Funció aov3 del package sasLM
library(sasLM)
aov3(anticos~dosis*vacuna,immuno_1)
## Response : anticos
##                 Df  Sum Sq Mean Sq F value    Pr(>F)    
## MODEL           11 155.052  14.096  3.3954  0.002813 ** 
##  dosis           2   8.284   4.142  0.9977  0.378963    
##  vacuna          3 142.953  47.651 11.4785 2.167e-05 ***
##  dosis:vacuna    6   5.375   0.896  0.2158  0.969257    
## RESIDUALS       35 145.297   4.151                      
## CORRECTED TOTAL 46 300.349                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aov3(anticos~vacuna*dosis,immuno_1)
## Response : anticos
##                 Df  Sum Sq Mean Sq F value    Pr(>F)    
## MODEL           11 155.052  14.096  3.3954  0.002813 ** 
##  vacuna          3 142.953  47.651 11.4785 2.167e-05 ***
##  dosis           2   8.284   4.142  0.9977  0.378963    
##  vacuna:dosis    6   5.375   0.896  0.2158  0.969257    
## RESIDUALS       35 145.297   4.151                      
## CORRECTED TOTAL 46 300.349                              
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.5.3.2 Amb el package car

## Funció Anova del package car
immuno_1.aov<-aov(anticos~vacuna*dosis, data=immuno_1,contrasts=list(vacuna=contr.sum,dosis=contr.sum))
car::Anova(immuno_1.aov,type=3)
Sum Sq Df F value Pr(>F)
(Intercept) 1170.163604 1 281.8765706 0.0000000
vacuna 142.952692 3 11.4784561 0.0000217
dosis 8.283744 2 0.9977210 0.3789631
vacuna:dosis 5.375372 6 0.2158091 0.9692567
Residuals 145.296667 35 NA NA
immuno_1.aov<-aov(anticos~dosis*vacuna, data=immuno_1,contrasts=list(vacuna=contr.sum,dosis=contr.sum))
car::Anova(immuno_1.aov,type=3)
Sum Sq Df F value Pr(>F)
(Intercept) 1170.163604 1 281.8765706 0.0000000
dosis 8.283744 2 0.9977210 0.3789631
vacuna 142.952692 3 11.4784561 0.0000217
dosis:vacuna 5.375372 6 0.2158091 0.9692567
Residuals 145.296667 35 NA NA

Pots trobar més informació en el següent enllaç

1.6 Dissenys amb més de dos factors

La teoria dels models lineals pot ser estesa a models amb més de dos factors, per exemple un model amb 3 factors fixos seria el següent

\[ Y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha \beta )_{ij}+ (\alpha \gamma )_{ik}+(\beta \gamma )_{jk}+(\alpha \beta \gamma)_{ijk}+\epsilon_{ijkl} \\ i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,c;\,\,\,\, l=1,...,n_{ijk} \]

amb

\[ \sum_{i=1}^a \alpha_i = 0 \,\, , \,\,\sum_{j=1}^b \beta_j = 0 \,\, , \,\,\sum_{k=1}^c \gamma_k = 0 \]

\[ \sum_{i=1}^a (\alpha \beta )_{ij} = 0 \,\, , \,\, \sum_{j=1}^b (\alpha \beta )_{ij} = 0 \,\, , \,\, \sum_{k=1}^c (\alpha \gamma )_{ik} = 0 \,\, , \mbox{ etc. } \]

\[ \sum_{i=1}^a \sum_{j=1}^b (\alpha \beta \gamma )_{ijk} = 0 \,\, , \,\, \mbox{ etc. } \]

Si dos o més dels factors son aleatoris per a alguns termes no hi haurà un quocient F adequat per contrastar la seva significació. El problema es torna quasi intractable si tenim múltiples factors aleatoris i un disseny no balancejat.

Com a mostra veiem la següent taula extreta del llibre de G. Quinn i M. Keough Experimental Design and Data Analysis for Biologists, Cambridge University Press, 2002.

A la taula \(\sigma^2_\alpha\) es refereix a la variància aportada si el factor A és aleatori o bé a \(\frac{\sum_{i=1}^a \alpha^2}{a-1}\) si el factor és fix.

1.6.1 Exemple i resolució amb R

colesterol<-read.table('../dades/cholesterol.csv', sep=';',dec=',',header=T,stringsAsFactors = T)
colesterol.aov<-aov(cholesterol~drug*sex*risk,data=colesterol)
anova(colesterol.aov)
Df Sum Sq Mean Sq F value Pr(>F)
drug 2 1.3952111 0.6976056 5.4661048 0.0065953
sex 1 0.0032000 0.0032000 0.0250737 0.8747154
risk 1 0.0000889 0.0000889 0.0006965 0.9790330
drug:sex 2 0.0079000 0.0039500 0.0309503 0.9695392
drug:risk 2 0.7256778 0.3628389 2.8430327 0.0661211
sex:risk 1 0.0060500 0.0060500 0.0474049 0.8283808
drug:sex:risk 2 0.0399000 0.0199500 0.1563187 0.8556338
Residuals 60 7.6574333 0.1276239 NA NA